The "Coulomb phase" in frustrated systems 



Christopher L. Henley 

Dept. of Physics, Cornell University 



Abstract. The "Coulomb phase" is an emergent state for lattice models 
(particularly highly frustrated antiferromagnets) which have local constraints that 
can be mapped to a divergence-free "flux" . The coarse-grained version of this flux 
or polarization behave analogously to electric or magnetic flelds; in particular, 
defects at which the local constraint is violated behave as effective charges with 
Coulomb interactions. I survey the derivation of the characteristic power-law 
correlation functions and the pinch-points in reciprocal space plots of diffuse 
scattering, as well as applications to magnetic relaxation, quantum-mechanical 
generalizations, phase transitions to long-range-ordered states, and the effects of 
disorder. 
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1. Introduction 

1.1. The basic idea 

A class of interesting lattice systems in solid state physics and statistical mechanics 
have ground states which (to a first approximation) are highly constrained yet highly 
degenerate - that is. more or less, the current definition [5] of a "highly frustrated" 
magnet. A fundamental problem in handling such systems theoretically is that of 
"navigation" among these states. Is there a way to label and enumerate the states so 
as to carry out a sum over them, and evaluate a partition function? Or if a particular 
state is somehow to be "selected" out of the ensemble, | how do we find that "needle 
in a haystack"? 

One of the standard answers is a coarse-graining, that discards most of the 
information in the configurations and keeps local averages of certain quantities, thereby 
converting the lattice problem into a continuum model. For this to be fruitful (and 
physically meaningful) the quantity being averaged ought to be conserved (It might 
also be the order parameter of some long range order, but our models are liquid- 
like in the first approximation - that is just paraphrasing the characterization in the 
first sentence.) In ordinary (off lattice) liquids, this is rather mundane: local particle 
densities and momentum densities. But in the lattice models I have in mind, we 
are back in the wilderness: the microscopic degrees of freedom (e.g. spins) usually 
don't have a conservation law, and (it being a lattice model) there is no momentum 
conservation. 

Fortunately, in quite a few cases, another kind of conservation is hidden in the 

first-order ensemble: a constraint on the total spin (or other degree of freedom) 
surrounding each lattice point. That allows us to map each microstate of local variables 
into a configuration of (weighted) arrows living on the bonds of the lattice, such that 
the signed sum of the arrow weights into every vertex (outwards minus inwards) is 
exactly zero in any allowed configuration. Such arrows are called "lattice fluxes" since 
this is exactly the zero-divergence condition an electric or magnetic flux would satisfy 
(in the absence of sources), if the field were constrained to lie along those lattice edges 
and its flux could only take discrete values. 

The desired emergent vector fleld P(r) is the coarse graining of this flux: [1, 2, 3, 4] 
that is, the mean value of the lattice fluxes over a volume centered at r much bigger 
than a lattice consstant, but much smaller than the system size; the divergence 
condition becomes V • P(r) = 0. Furthermore it turns out [Sec. 3] that the effective 
coarse-grained free energy has the form J d'^r^\P\^, exactly the form of the field 
energy of an electric (or magnetic) field. 

There are many ways to apply this analogy to find the long-distance behaviors 
of oTir constrained model. Since the probabilities according to this free energy are 
Gaussian, one can compute practically any desired expectation. In particular, the 
"spin" correlations depend on separation R with the functional form of a dipole- 
dipole interaction, proportional to in d dimensions. It was somewhat siuprising 

to find such a slow decay (with a divergent correlation length) in such a liquid-like, 
maximally random system. Such states (in three dimensions) have acquired the name 
"Coulomb phase" . 

The conditions for a Coulomb phase are that 

I For example, the classical state which optimizes the energy gain from some perturbation term in 
a classical Hamiltonian, or from quantum fluctuations in a quantum Hamiltonian. 
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(CI) each variable can be mapped to a signed fliix Pj running along bond i; 

(C2) the variables obey hard constraints, such that the sum of the (incoming) fluxes 

at each (parent) vertex is zero; 

(C3) the system is in a highly disordered phase, without any long range ordered pattern. 
This may be called "liquid-like" to express the disorder coupled with strong local 
correlations [implicit in (C2)]. 

1.2. Rearrangements 

What local rearrangements are permitted by the flux constraint? If I flip the variable 
on a bond adjacent to parent lattice site a, so as to change its flux from incoming to 
outgoing, I must flip the flux in the opposite way on one of the other bonds. Thus, 
the natural rearrangements (either in a simulation, or a real system) are entire loops, 
which are sometimes called "Dirac strings" . As one traverses the string in a particular 
direction, the sense of the flux arrows is always the same with respect to the walking 
direction; the local rearrangment reverses this sense, from always forwards to always 
backwards or vice versa. Thus, if the string extends across the whole system, the 
rearrangement changes P, but if the string closes, then P is unchanged. [An example 
of such an update is shown later in Fig. 4(a).] 

1.3. Outline of the article 

The aim of this review is to survey the lattices and models (and real materials) in 
which a Coulomb phase are found (Sec. 2); to walk through the derivation of its power- 
law correlations (Sec. 3); and to highlight some interesting ways that Coulomb-phase 
ideas have been deployed to solve problems in frustrated systems. These include, 
topological defects that (in the case of dipolar spin ice) are (emergent) magnetic 
monopoles (Sec. 4); dynamics (Sec. 5); quantum- mechanical generalizations (Sec. 6); 
transitions to ordered phases (Sec. 7); and quenched disorder (Sec. 8). 

Although this review touches on many aspects of frustrated spin models (gauge 
theories, disorder, quantum dimer models, dynamics), it does not aim or claim to 
review any of these major topics, except for particular instances that happen to be 
tractable using the Coulomb-phase notions. Indeed, even concerning the Coulomb- 
phase developments to date, I have not tried to exhaustively survey all papers - I 
only try to represent each aspect of the topic in the sections. My choice of particular 
results to highlight is driven by the motives of presentation, and does not always imply 
priority or importance. 

2. Examples 

How can wc realize a Coulomb phase? This breaks up into two questions: which 
models, on which lattices, have the requisite constraints? (Sec. 2.1). And, which 
physical systems realize such models (the rest of this section). Of course, though two 
physical systems may realize mathematically equivalent models, quite likely different 
kinds of quantities are experimentally accessible in the respective systems. 

I will also point out that one might view the constrained system as either a T = 
or T = oo limit of some statistical mechanics model. If we consider an enlarged 
configuration space in which the constraint is not forced, but its violation costs energy 
(as I do in Sec. 4), the constrained model is the T = limit. On the other hand, within 
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the constrained ensemble we might consider adding a Hamiltonian which breaks the 

degeneracy among states examples arc kagome ice (Sec. 2.5.1) and the perturbations 
driving transitions in Sec. 7 - in which case the basic ensemble is the T = oo limit. 

2.1. Lattices 

What models give rise to Coulomb phases? First I will talk about the lattices, and 
then the degrees of freedom. 

2.1.1. Parent and medial lattices First of all, we always have a "parent" lattice B 
which is bipartite, i.e. one that can be partitioned into (equivalent) sublatticcs of even 
and odd nodes such that every bond connects an even node to an odd node. Degrees 
of freedom - call them "variables" for short - live on these bonds, i.e. on sites of 
"medial" lattice C consisting of the bond midpoints. (I will call sites of the parent 
lattice "nodes" simply to help keep them straight from "sites" of the medial lattice.) 

Conversely, I call B the "premedial" lattice of jC. (In the literature, B is often 
sloppily called the "dual" lattice, which properly means something different. For 
example, the premedial lattice of the kagome lattice is the honeycomb, but the dual 
of the kagome is the dice lattice.) 

See Tabic 1. The most important three-dimensional example is the pyrochlore 
lattice (diamond as parent lattice). The simplest realization is the "i?" sublattice 
(octahedral sites) of the spinel (usually oxide) structure. The other important 
realization is the pyrochlore crystal structure. (That is a large family of oxides, which 
contain two interpenetrating pyrochlore lattices, each occupied by a different species 
of cation.) 

An example which deserves more attention is the "half-garnet" lattice. The 
magnetic lattice in a garnet consists of two interpenetrating copies of this lattice. 
§ The parent lattice is the "Laves graph of degree three", in which each vertex has 
three bonds forming 120° angles, and its Bravais lattice is bcc; thus, its symmetry is 
just as high as the pyrochlore's. 

I also included the "octahedral" lattice, which gets reinvented from time to time 
as a toy model because its parent lattice is simple cubic, suitable for simple-minded 
theorists. Finally, the "sandwich" lattice consists of two kagome layers linked by an 
additional triangular layer, and models the antiferromagnet SrCr8_2;Ga4+a;Oi9), so it 
is essentially two dimensional. There are experimental realizations for most of these 
lattices. 

For pedagogical purposes, I will use two-dimensional lattices in all figures, even 
though (see Sec. 3.3) these are not quite bona fide Coulomb phases. The reader should 
view them somewhat more as analogies (to d — 'S) rather than as examples. 

2.1.2. Degrees of freedom The most common models reduce to either ice- models or 
dimer models. An ice model is defined on a parent lattice with even coordination 
number (also called "six- vertex" model in the usual case Zs = 4.) Every edge 
carries an arrow, and at every vertex the "ice rule" constraint is obeyed, namely half 
the vertices point in and half point out [14, 15]. This arrow obviously is the flux. [See 
Fig. 1(a).] 

§ The "hyperkagome" lattice [87, 88] is equivalent to the half-garnet lattice, so long as only first 
neighbor bonds are talsen into account, but has less symmetry: one primitive cell of hyperkagome is 
two primitive cells of half-garnet. 
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Table 1. Parent lattices and medial lattices. The table shows the spatial 
dimension d, the parent lattice with its coordination Zg, and Bravais lattice, 
the medial lattice with its coordination Z^, and the number of sites per cell. The 
length ^loop of the shortest loop (same in either lattice) is also given. Finally, a 
reference for the lattice is given. The clusters formed around each node of the 
parent lattice have sites (forming triangles, tetrahedra, or octahedra). 



d 


parent 


Zb 


Brav. 


medial 


Zc 


sites 


-^loop 


Rcf. 




lattice B 




latt. 
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Laves graph 
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b.c.c. 


half-garnet 


4 


4? 


10 


[13] 



A dimer means an object which covers two nodes of the parent lattice. Unless 
otherwise specified, a dimer covering also satisfies the condition that every node is 
covered by exactly one dimer: dimers never overlap, and no node is left uncovered. 
To define the flux, we draw an arrow of weight Zg — 1 along each dimer-occupied 
edge, pointing from the even to the odd node; and arrows of weight 1 along each 
unoccupied edges, pointing in the opposite direction. (Obviously a different overall 
normalization of the flux could be used.) Illustrations of this weighting may be found 
in Fig. Al(a,b,d). 



(a). (b). (c). (d). 




Figure 1. Mappings of the ice model on the diamond lattice (green edges), (a). 
Polarization arrows of the abstract ice model, (b). Water ice (c). Compound 
of species A and B (d). Ising ground state. Each configuration shown in (b,c,d) 
maps to the arrow pattern in (a). 



2.2. Water ices 

I now turn to the different kinds of physical realization, starting with real ice. In 
crystalline II2O, the O atoms form a diamond lattice with hydrogen bonds to all four 
neighbors. || Two of the these adjacent H's are covalently bonded to the O in question 
(as H2O), while the other two H's belong to the neighboring O's. [Fig. 1(b)]. The 
variables representing a configuration of the H's are arrows along the lattice bonds, and 

II The usual form of water ice uses the hexagonal diamond lattice. The cubic diamond lattice used 
in the model has very similar behaviors, but is nicer theoretically owing to its higher symmetry. 
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clearly satisfy the "ice rule" . This is the first highly frustrated model to be analyzed, 
by Bernal [14] and Pauling [15] in the 1930s. 

There are many structures (often ferroelectrics) having hydrogen bonds such that 
the low energy structures realize some kind of ice rules. 

2.3. Lattice- gas orders 

An important way to realize constraints is when some kind of mutually repelling 
particle occupies the medial lattice sites. If we constrain the overall filling n by these 
particles to certain rational fractions, and if the interactions are only nearest neighbor, 
it is easy to see the ground states are those in which every triangle or tetrahedron 
has the same filling. If have two species A and B which attract each other more than 
their own kind, with filling n of species A and 1 — n of species B, as suggested for 
CsNiCrFe [16, 17] we get exactly the same ensemble. 

When n = 1/2, we can map the configurations to an ice model [Fig. 1(c)]: each 
occupied site becomes an even-to-odd arrow on the corresponding bond of the parent 
lattice; each vacant site becomes a bond in the opposite direction [18]. Another special 
filling is n = 1/Zg, in which case we map the occupied sites to dimers. 

One realization is in oxides where the cations have a mix of two valence states, 
corresponding to species A and B; an example of this is magnetite Fe304, which has 
an equal mixture of Fe+^ and Fe+^ on the spinel-B-sites [18]. 

2.3.1. Heavy- electron spinel LiV20i Another example is the metallic spinel 
LiV204 [19], which exhibits hcavy-fermion behaviors despite the absence of / 
electrons. With ion charges Li+1 and 0~^), charge balance demands (on the vanadium 
occupied B sites) n = 1/2 of and the rest of V"'"''. Fulde and collaborators [19, 20] 
proposed quantum-mechanical models in which the correlated electrons can hop with 
amplitude t on the frustrated lattice. ^ They work from the limit of strong Coulomb 
repulsion, modeled discretely in the spirit of the Hubbard model; an on-site term 
effectively constrains the occupancy to at most one electron per site, plus the inter- 
site repulsion V which is responsible for the flux constraint. 

2.4. Antiferromagnets 

The remaining kinds of realization are magnetic. 

2.4.1. Relation of antiferromagnetic Hamiltonian to constraint Define to be the 
total spin on the cluster surrounding site a of the parent lattice: 




(2.1) 



(denoted by "i e a"). We want the Hamiltonian to constrain 
La =0 

for all a. Next, say the Hamiltonian can be written as 



(2.2) 




(2.3) 



a 



^ Since the orbitals are actually d orbitals, I have glossed over several complications of orbital 
degeneracy and orbital dependent hopping amplitude, as well as spin. 
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Manifestly, this gives what we wanted: the classical ground states are any configuration 
satisfying (2.2), - the net spin in every cluster is zero. Furthermore, all such 
configurations are degenerate. 

But in fact, if one simply expands the square in (2.3), one gets an antiferromagnet 
with coupling Jy — Jgpin if i and j belong to the same cluster, and zero otherwise. In 
the example lattices (all except the octahedral) where the every cluster is a triangle or 
tetrahedron, that just means the nearest neighbors: the ground states of the nearest- 
neighbor antiferromagnet are highly degenerate, but every state satisfies the constraint 
(2.2) around every a. 

The ground state ensemble of a system obeying constraints and lacking any kind 
of long range order may in general be called a "cooperative paramagnct" (sometimes 
called "classical spin liquid" which means the same thing.) It bears the same relation 
to the usual high temperature paramagnet 

Next I explain the three different situations in which antiferromagnets can form 
Coulomb phases [3]. 

2.4-2. Ising antiferromagnet The ground states of the pyrochlore lattice Ising 
antiferromagnet have two up and two down spins in every tetrahedron. The Ising 
model is not realistic in its own right for three-dimensional lattices (mainly because 
their symmetry is incompatible with the special axis of the Ising spins), but many other 
models map to it. For example, the pyrochlore Ising antiferromagnet 's ground states 
map 1-to-l onto those of diamond-lattice ice model [18, 21]. Each spin ti = +1(— 1) 
maps to an arrow pointing along the corresponding diamond lattice edge, in the 
positive (negative) sense from the even to the odd vertex. [See Fig. 1(d).] 

Notice that, if we turn on an external field tuned to the appropriate size, the 
ground states (on any of our lattices) have — 1 spins up and one down on every 
cluster, so this maps to a dimer covering. 

2.4-3. Isotropic Heisenberg model: "classical spin liquid" This state arises when the 
spins arc classical vectors (or can be treated as such), ?/ all states in the contimiOTis 
manifold satisfying (2.2) are more or less equally likely. In this phase, every vector 
component of the spins satisfies the flux constraint; thus, the polarization field has 
indices not only for directions in space, but also for the three spin components. 

2.4- 4- Isotropic Heisenberg model: "Emergent discrete spins" This state can arise 
in the classical models (mainly the triangle-based lattices) that manifest "order by 
disorder" in the sense of Moessner and Chalker [6, 7]; it also arises in all quantum 
models with sufficiently large spin length S, when T is low enough that spin-wave 
energies favoring coUinear or coplanar states are important [22, 23, 24]. We get a 
discrete subset of spin states in which all spins are coUinear along a given axis (or 
coplanar when the clusters are triangles). The ensemble of this discrete subset is some 
sort of Ising model or coloring model, still satisfying the constraint (2.2). 

2.5- Spin ice 

In "spin ice" pyrochlore magnets Dy2Ti207 andHo2Ti207 local (111) spin anisotropics 
(additional to (2.3)) reduce the ground state manifold to effective Ising states [25, 26], 
as reviewed in [27]. The easy axis is the bond direction Uj: Si = fjUj with ti = ±1. 
The actual interactions are ferromagnetic, Jp < 0- Since Um ■ Um' = for 
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neighboring spins, the Hamiltonian in terms of {ti} reduces to an Ising model with 
antiferromagnetic J = —Jp/3. In terms of the original model, the groimd states of 
this ensemble literally implement the ice rules, i.e. in each tetrahedron (surrounding 
a parent lattice node), two spins point in and two point out. 

2.5.1. Kagame ice When a (not too strong) magnetic field is placed on spin ice 
along (say) the [111] direction, it selects out a subset (still with extensive entropy) of 
ground states. Namely, every spin on a bond of the parent lattice in the 111 direction 
is parallel to the field. The three remaining sublattices are still free to fluctuate; but 
they form disconnected kagome layers, so we now have a stack of independent two 
dimensional systems, hence this system (realized experimentally) is called "kagome 
ice" [28]. The parent lattice of each layer is a honeycomb lattice; the spin constraint 
in the layer is two in/one out on (say) the even nodes, and oppositely on the odd nodes. 
These configurations map exactly to those of a dimer model on the honeycomb lattice. 
(The dimer positions correspond to the single spin with its in/out sense opposite to 
the others in its triangle.) 

2.5.2. Dipolar spin ice The actual spin ice materials Ho2Ti207 and Dy2Ti207 are 
well approximated as having nothing but (long-ranged) dipolar spin interactions, 
rather than nearest-neighbor ones. Although this model is clearly related to the 
"Coulomb phase" , I feel it is largely an independent paradigm with its own concepts 
that are different from the (entropic) Coulomb phase that most of this review is about. 

The key fact is that all microstates satisfying the ice rules have very nearly 
the same energy [29]; this was first appreciated numerically [30] (after the technical 
problems of simulating long-range interactions in periodic boxes were sorted out). The 
simple explanation was given in [31]. It's a good approximation - for the far field, 
anyhow - to replace each point dipole by a pair (so-called "dumbbell") of opposite 
effective charges ±q having the same dipole moment jj,, and we might as well separate 
them by the bond length d of the parent lattice, thus fi = qd. So, all the effective 
charges sit on parent lattice nodes, and each node gets contributions from four dipoles 
(see Fig. 2). Indeed, so long as the spin configuration satisfies the ice rules, those 
four contributions add to zero on every node, and (within this approximation) the 
interaction energy is zero (I omitted a configuration-independent constant representing 
the interactions between the two charges forming each dipole). In principle, small 
multipole interactions, representing the difference between the actual interactions of 
point dipoles and those of the effective charges, will break the degeneracy and lead to 
ordering at much lower temperatures. 

3. Coarse-graining, Fourier mode fluctuations, and long-range correlations 

In this section, I present the steps leading to power-law correlations using the 
framework of a continuum theory, in the ideas appear more transparently. 

3.1. Polarization field and effective free energy 

The polarization field is the key object of the Coulomb phase. As already laid out, our 
system has an extensive entropy of ground states (degenerate at this order) . How to 
handle this? Recall the discrete fluxes Pj, which we defined along every bond of the 
parent lattice. Define a (whole-system) polarization density P = J2i Pj/volume. This 
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(a). (b). 




Figure 2. Dipolar spin ice and monopoles. (a). A configuration obeying ice 
rules, except for defects of charge Q = +2 and —2 (marked), (b). When each 
dipole is replaced by a pair of charges (red + and — symbols) at the ends of its 
bond, these cancel out at every node of the parent lattice where the ice rule is 
obeyed (open blue circles), but add to a total proportional to Q at defect nodes 
(shaded purple circles). 



is a good choice for writing thermodynamic functions, since the ground-state entropy 
density evidently depends on it. 

Recall the point of Sec. 1.2, that rearrangements (within the flux constraints) 
must occur along strings which follow flux arrows. Now, if P is large, most fluxes have 
a positive component in that direction and the typical string crosses the system - a 
closed string would require half of its fluxes to have a backwards component, but such 
a local gathering of reversed fluxes is unlikely. On the other hand, when P is near zero, 
there are many closed strings in any configuration, and hence many rearrangements 
which preserve P. It's clear, then, that the number of configurations A/'(P) with a 
given polarization is maximum at P = and goes to zero as P approaches its saturated 
value. So let's define an entropy density 

,(P) . li„ (3.1) 

which is maximum at zero. 

Next, so long as the ensemble is a "liquid" lacking long-range orders, the local 
polarization (if we divide the system into smaller boxes) is fluctuating, and not too 
strongly correlated from box to box. Then the Central Limit Theorem says the 
ensemble probability (= number of ground states for a given P) in a system of large 
volume V approaches a Gaussian form 

AA(P) oc exp(-|P|V2(7p2) (3.2) 

with a variance dp = 1/KV for some K. Comparing with (3.1), we see 

s(P) « so - |PP (3.3) 

for small P. (All this was a completely standard argument from basic statistical 
mechanics; it is the same reason the free energy is proportional to (magnetization)'^ in 
a paramagnet, or that one assumes analyticity in the Landau free energy functional.) 

Next let's consider the spatial fluctuations of P, which are more interesting (and 
more measurable!) - than the functional form of s(P). This requires defining a 
spatially varying polarization field P(r); it is a coarse-graining, i.e. the average of 
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discrete polarizations over some neighborhood of r (one much larger than the lattice 

constant but much smaller than the system size.) Wc assume P(r) varies smoothly. 
Corresponding to the discrete flux constraint, P(r) satisfies a divergence constraint 

V • P(r) = (3.4) 

like a magnetic field without monopoles. 

The total free energy (arising entirely from entropy) is the sum of those in the 
boxes or averaging volumes into which we can divide the system, hence 

Ftot({P(r)})/T = const + j d''r^K\P{r)\\ (3.5) 

Eqs. (3.5) and (3.4) look, respectively, like the field energy of a magnetic (or 

electric) field, and its divergence constraint, in the absence of monopoles (or charges). 
That electrostatic (or magnetostatic) analogy is fruitful, and is why this state was 
dubbed "Coulomb phase." 

In the case of dipolar spin ice, since the fluxes are parallel to real moments, the 
polarization is proportional to the real magnetization M. 

M = ^P. (3.6) 

Furthermore, if a region has a net polarization, we have a field energy of form 
(3.5), but now K = Ho^i^, where fj,a is the permeability of free space, in the hmit 
T = 0: it is purely energetic rather than entropic. At T > 0, the cntropic elasticity 
gives an correction to the permeability. (In water ice, which has long range electric 
dipole interactions, the analogous contribution to the dielectric constant has long been 
known [32].) 



3.2. Pseudodipolar correlations and structure factor 

The standard way to evaluate correlations is to transform to Fourier space. Eq. (3.4) 
gives 

q-P(q) = (3.7) 
so a naive use of equipartition would give (P^(— q)P!y(q)) = {1/K)5fj,v. so (3.5) gives 

q 

where P_L(q) refers to the components of P(q) satisfying (3.7). The fluctuations are 
gotten by first writing the (trivial) result of equipartition for an unrestricted P(q), 
and then projecting to obtain the transverse part: 

S,M) ^ (P^(-q)P.(q')) = - (3.9) 

3.2.1. Diffraction consequences A physical observable <I'(r) usually has a contribution 
proportional to P, but possibly the correspondence is modulated e.g. by alternating 
signs, so we write 

$„(r) = (...) + ^c„^(r)P^(r). (3.10) 
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where {ca^(r)} is a matrix of coefficients, with the symmetry of the lattice, and indices 

ab refer to possible i components of the observable. It follows that the observable 
structure factor, measured in diffraction, behaves as + 



where /a6^^(Q) is some sort of structure factor derived from the Cafi^s, and Q are 
reciprocal lattice vectors. 

Two important features may be noted. First, the polarization constraint says 
that a certain linear combination of fluxes is exactly zero (their projection on the 
inward direction at a node of the parent lattice). That translates into a strict 
zero of the appropriate structure factor, at a certain wavevector. In the case of an 
antiferromagnet, or spin ice, the constraint is that magnetization is exactly zero in 
each cluster. Hence, we find a very low scattering near the zone center (the diffuse 
intensity grows as a high power of |q|); also, the intensity has zeroes along every 
direction that is a symmetry axis of all clusters. 

The second and more striking feature comes from the second term in (3.9). 
Although not divergent, it is singular: the ratio has a different limit at q = 0, 
depending on the direction of approach. In reciprocal space, this has the characteristic 
shape of a "pinch point", at which the contours of equal intensity have a roughly 
triangular shape. The form factor in (3.11) translates these to reciprocal lattice vectors 
Q other than zero, as seen in date from [33], shown in Fig. 3. (Earlier experiments [34] 
gave less clearcut images: polarized neutron diffraction is needed to separate the 
contributions from different spin components.) 

I discussed in [3] the way that other kinds of hard constraint also produce sharp 
features in reciprocal space. For example, if you constrain a lattice gas on an fee 
lattice so that every particle has an equal number of occupied neighbors, you get 
rings [35]. (Similar features are observed in quasicrystals and ascribed to local tiling 
constraints [36]). This constraint superficially resembles that of a coulomb phase as 
laid out in Sec. 2.3; the outcome is different because the fee lattice isn't the medial 
lattice of a bipartite lattice.) 

3.2.2. Real space correlations Fourier transforming (3.9) back to direct space gives 



at large separations r in d dimensions, where f = r/)r) and C3 = 47r. So the 
correlations, which one naively expected to be exponentially decaying in this "liquid- 
like" state, are instead power-law-decaying i.e. criticalr\SkB, They have the spatial 
dependence of a dipole-dipole interaction. 

This criticality was appreciated in the ice model as early as 1973, being detected 
originally in a simulation [37]. The universal explanation (above) how dipolar 
correlations arise from (3.5) with the divergence condition, was first put forward in 
1981 to explain experiments on two-dimensional ice-like systems [1]. It was noticed in 
the context of antiferromagnets in [38] and (most infiuentially) announced as a general 
idea (initially for dimer coverings) by Huse et al [2]. 

We can apply these ideas to the nonlocal spin susceptibility near a defect. 
Consider a vector-spin antiferromagnet on one of the lattices in Table 1, where the 

+ If "(...)" in Eq. (3.10) is a function of local fluctuating variables that has no conservation law, then 
"(...)" in (3.11) adds a non-singular diffuse contribution on top of the features of <Snii/(q). 



<S* (q) ^ (I.„(q)$,(-q)) = (...) + ^ /a(,;.^(Q) V(q - Q). (3.11) 




(3.12) 
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spins are diluted (this corresponds to bond dilution on the parent lattice). At a 
dangling node, one that has just one neighbor, the flux constraint cannot possibly be 
satisfied. Such "orphan spins" respond to external fields like free moments, and in line 
with Eqs. (3.12) [or (4.1), below] their perturbation of the surrounding spins decays 
as a power law (with oscillations depending on how spins map to polarizations) . This 
can be modeled within a Coulomb-phase framework, [39] and used to explain NMR 
observations in SrCr9pGai2-9pOi9. 

3.2.3. Calculation for specific lattices So far we only obtained the asymptotic 
functional form near "pinch points" : more practically, one would like to model 
(approximately) the diffuse scattering over the entire Brillouin zone for a specific 
model. For the spin models of Sec. 2.4, a systematic way is to allow n ^ oo where 
n is the number of spin components. This "large-n" approach is a standard trick 
of statistical mechanics, since (in that limit) the constraint of unit length becomes 
irrelevant and we have a linear problem [4, 40, 41]. This turns out to be an good 
approximation even for n = 1. 

An ad- hoc "maximum- likelihood" approach, [3], which can be applied to any 
model of our class, is to consider the fluxes Pi (representing the model variables) to be 
rea/ numbers, subject to two kinds of constraint: (i) that take certain discrete values 
(ii) the usual divergence constraint. We then replace (i) by a weighting exp(— iPj^/cr^), 
where cr^ is chosen to give the correct variance when the {Pi\ are chosen at random 
(and unconstrained by the divergence condition) . The cluster- variational approach [42] 
is related. 

Finally, Villain predicted correlations in ice using a clever random-walk 
approximation to a series expansion [43], which has not yet been reconsidered in 
the literature. I would suspect this amounts to a Bethe-lattice approximation (i.e. 
neglecting the existence of loops in the lattice). 




-2-1012 -2-1012 



Figures. Pinch point in diffraction. From [33]. Spin-flip component of polarized 
neutron diffuse scattering from the spin-ice pyrochlore compound Ho2Ti207, in 
the (hhl) plane of reciprocal space. Left side, experimental result; right side, 
Monte Carlo simulation. The added arrows highlight the "pinch point" features 
at the wavevectors labeled. 
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3.3. Height models: two-dimensional Coulomb phases 

It is even easier to realize a flux constraint in d = 2 than in = 3, and such models 

exhibit all the phenomena mentioned hero. However, they exhibit additional, even 
more striking, behaviors which are peculiar to two dimensions, and for that reason 
they fall outside the scope of this article. Here I will just summarize the differences; 
such "height models" call for a review paper of their own [44, 45, 46, 47, 48, 49]. 

To see what's different, recall that a divergence- free field such as our P(r) can 
always be written as a curl of a vector potential, P(r) = V x A(r) . In three dimensions, 
A is ill defined due to its gauge freedom; but in two dimensions, ^(r) has only 
one component (normal to the plane) and is an ordinary potential, uniquely defined 
modulo an additive constant. [To see this, note that Prot = {—Py,Px) carries the same 
information as P but satisfies V x Pj-ot = 0, hence Prot = defines a potential.] 

We can think of this as mapping a configuration of the variables to one of a 
crystal interface [44, 45] with profile z = h{x,y) = A{x,y). The "Coulomb phase" 
corresponds to the "rough" phase of such an interface model with the well-known 
form i^tot = / d'^r^K\Vh{r)\'^ . However, because h{r) is uniquely defined, a generic 
physical variable (written in terms of field variables) is not only a linear combination 
of Vh terms, but also has periodic terms of form cos(27rm/i/aj^). [The repeat offset 
a± depends on the model and the definition of h{r).] It can be shown that: 

(i) This leads, in the rough phase, to critical correlations with a parameter (and 
temperature) dependent exponent oc T/K. 

(ii) An un-rougliening transition occurs when T/K decreases past a universal ratio, 
in which the system locks to a particular value of A, corresponding to long-range 
order of the model variables; 

(iii) Thermal excited "defect charges" (see Sec. 4) do not necessarily destroy the 

critical phase; rather, they stay bound in pairs for T/K less than a universal 
critical ratio; the unbinding is the "Kosterlitz-Thouless" transition familiar in 
the XY model. 

The critical correlations (i) imply singularities in reciprocal space, generally 
displaced at different positions in the Brillouin zone than the pinch points. Such 
singularities (in the height field) were called "zone-boundary singularities" in Refs. [47, 
48, 49] or (in the structure factor) "pi-ons" in Ref. [50]. 

4. Pseudo-charge defects in Coulomb phases 

If we change our rules to allow configurations that violate the flux constraint, merely 
including a Hamiltonian which strongly penalizes the violations so they are dilute at 
low temperatures - what do we get? Such a defect may be labeled by its "charge" 
Q, equal to the net (nonzero) flux in the outward sense (at a parent-lattice vertex). 
This is a "charge" in the sense of Gauss's law, as the net flux through any surrounding 
surface must equal Q. That means the "charge" is conserved in time: the net charge in 
some volume can't be changed except by moving it across the boundary, and defects 
can only be created in pairs of opposite charge. Their detectability from distant 
measurements, as well as the conservation properties, are characteristic of topological 
defects. 
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More interesting, in a Coulomb phase, the effective potential between defects 
must have the form of Coulomb's law: 

Fi„t(ri,r2)/T=-^^i^ (4.1) 

(in d = 3). The interaction is defined by integrating the partition function, conditional 
that defect charges Qi and Q2 are placed at ri and r2, and by requiring the result to 
be proportional to exp{—Fint/T). That is essentially how one obtains the Coulomb 
potential from the field energy in electro- or magnetostatics, and we must get the same 
result here since the field free energy has exactly the same form. In d = 3 (but not 
necessarily in d = 2), our "charges" are deconfined i.e. a defect/anti-defect pair (in a 
large system) will separate and have independent positions. By contrast, if these same 
defects are present in a phase with long-range order of the variables, pulling apart the 
defects now carries a cost proportional to |ri — r2|: in that case, defects are "confined" 
the same way that quarks are. 

Some authors have emphasized the notion of "Dirac string" meaning the trail 
of fluxes that got flipped as you pulled apart a defect/antidefect pair. It should 
be realized that this is a nebulous and not very helpful notion when applied in the 
Coulomb phase proper (with smallish polarization), since the string's path is not well 
defined. That is, there are many different, equally good ways to represent a defect 
pair by (i) taking some configuration from the undefected ensemble and (ii) flipping 
the fluxes along a string. It is only in an ordered phase (see preceding paragraph), or 
near the limit of maximum polarization, that the Dirac string has a clear meaning. 

When the model is a lattice gas of neighbor-repelling particles, that map to a 
dimcr covering on the parent lattice, a simple way to create a defect/antidefect pair 
is to remove one particle or dimer (the total flux at the parent sites that the dimer 
spanned gets increased/decreased by the flux associated with a dimer). Evidently 
each of the two defects is carrying an effective particle number of 1/2, making this an 
elementary example of "fractionalization" [51]. 

In actual water ice, OH~and H3O+ ions (if the hydrogen bond network is 
unbroken) are topological defects of charge Q = ±2 like the ones we were describing 
- except they aren't very mobile, since proton transfer is slow at icy temperatures. A 
distinct kind of defect, the Bjerrum defect, is formed when there are no protons, or 
two, along a given bond (maps to a nonmagnetic or doubly magnetic impurity in a 
spin system). The Bjerrum defects have charge Q = ±1 and are mobile. 

4-1- Thermal consequences 

If the creation energy is Eq for a defect of charge Q, the defect density will 
behave as uq oc exp{~EQ/T). The specific heat is affected proportionately. [As 
in semiconductors, the exponent is Eq/T even though each defect pair costs 2Eq, due 
to entropy , as the defect locations are independent.] 

The defect "charges" are in the same situation as ions in a plasma (or in water), 
or as carriers in a compensated semiconductor: we get Debye screening of the fields. 
In consequence, the "Coulomb" effective interaction, (4.1), acquires a decaying factor 
exp(— K;|ri — r2|), where 1/k is the Debye screening length, given by 



=2^. (4.2) 
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The pseudodipolar polarization correlations [Eq. (3.12) acquire the same exponential 
screening factor. In Fourier space, this corresponds to replacing |qp |qp + in 
the denominator of Eq. (3.9). 

How thermal excitation plays out in real systems depends on details of the model. 
In "spin ice", the spin configurations are not inherently constrained to obey the 
constraint, so at T > there is indeed a thermal concentration of defects. The 
cost 

Eq = JQ^ (4.3) 

is 2 J for the minimum charge Q = ±2. Following (4.2), then, we expect a 
K oc exp(— J/T), and this agrees with the experimental fit to the structure factor [33]. 
Around the reciprocal-space point where the defect-free diffuse scattering has a zero 
and flat basin, due to the exact cancellation in each tetrahedron, we expect an increase 
of diffuse scattering cx nq cx exp(— 2J/T) as is also observed [33]. 

On the other hand, say one has a dimer model on a triangular lattice and that 
dimer constraint is absolute. It still is not a Coulomb phase [51] (since the parent 
lattice isn't bipartite). However if a very strong external strain field biases the dimer 
orientations such that one orientation is excluded, the remaining (rhombic) lattice 
is bipartite and develops a Coulomb phase. (Exactly the same thing happens with 
dimers on the fee lattice. 

The experimental example of "kagome ice" , where the external bias is the 
magnetic field, is similar: in that case, the bias causes a two-dimensional Coulomb 
phase to emerge from a three-dimensional one.) 

Now, if the bias is reduced (but still strong) so that an occasional "wrong" 
orientation appears by thermal activation, that merely creates a dipole in the rhombic 
lattice's polarization field. As we also know from electrostatics, an ensemble of dipoles 
changes the dielectric constant (or its analog 1/K,m our case) but does not give Debye 
screening, so we continue to have power law correlations. 

4-2. "Magnetic monopoles" in dipolar spin ice 

In the case of dipolar spin ice (recall Sec. 2.5.2), what happens at a node where the 
ice rules are violated, so the net flux there is Ql The same construction outlined in 
Sec. 2.5.2, whereby point dipoles get replaced by pairs of effective charges ±q, must 
put a net (magnetic) charge of Q on every such defect node, just as it puts zero charge 
on every other node. (See Fig. 2). The total energy of such a state is now given 
by the Coulomb interactions of these defect charges: they are emergent magnetic 
monopoles [31]. 

We should be clear in what sense this is a monopole. Microscopically, of course, 
the laws of nature continue to rule out monopoles. Thus, near to any defect site, 
in between the bond directions where there is an excess of (say) outward flux, there 
must be other directions where there is a "counterflow" of inward flux. This seems 
rather analogous to the current associated with the Bogoliubov quasiparticle in a 
BCS superconductor. Its motion across the sample indeed causes a charge transport, 
proportional to the difference between the fractions of electron and hole making up 
the quasiparticle. Yet, while traveling in the bulk, the quasiparticle actually can have 
no charge [52] (in consequence, there, of the Meissner effect). 

Notice that although the "charge" interactions in dipolar spin ice look formally 
like those in entropic Coulomb phases, the physics is quite different. In the entropic 
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case, the energy cost is only the core costs of making the defects; the effective Coulomb 
interaction is wholly entropic, and emerges only as wc allow the system in equilibrium 
to explore all those many microstates. By contrast, in the dipolar case, every one of 
the (many) microstates having defects of +Q and —Q at ri and r2 has very nearly 
the same spin interaction energy, dependent only on |ri — r2|. 

A way to note the difference is that for an ice model in d = 2, defect charges or 
monopoles separated by R feel a In i? potential in the entropic case (just consider the 
Gauss's law for the model flux, confined to the plane) whereas in the dipolar case they 
have the usual 1/R interaction characteristic of three dimensions. Another situation 
that sharply illustrates the difference is "kagome ice" (see Sec. 2.5.1). So long as the 
ice rules arc satisfied, that system breaks up into decoupled two-dimensional layers. 
But when charge defects are introduced in the dipolar system, their interaction is 
isotropic - it is the same whether the two defects are in the same or in different layers 
(even though each defect is still free only to move within a layer). By contrast, in 
the kagome ice phase of model spin ice with nearest-neighbor couplings, an interlayer 
defect interaction exists at all only to the extent that spins linking the layers have 
some fluctuations and mediate it. 

One might reckon that the short-range models, in which (entropic) Coulomb 
behavior emerges from a collective state, is less trivial than the dipolar model, in 
which we have a bare magnetic energy. After all, our defect monopoles arc not terribly 
different from the emergent monopoles found at the end of a long thin bar magnet, 
as we were taught in introductory magnetostatics. The important difference, and the 
reason for all the excitement is that - unlike the bar magnet poles - 

(i) our monopoles move freely in response to forces; 

(ii) their magnetic charges are quantized (the value depends on material properties), 
the defect monopoles are mobile 

(iii) A variety of cute experiments is possible which depend on the monopoles. 

I would reserve the term "observation" of a monopole for an experiment isolating 
a single one. One might observe quantized jumps in the induced current around a 
conducting ring embedded in a sample whenever a monopole passes through it [31]. 
Alternatively, a monopole just inside the sample's surface creates a characteristic 1/R 
field outside it, which might be detectable by a scanning magnetic force microscope. 

The experiments so far are instead thermodynamic and transport measurements, 
the interpretation of which depends on monopoles being the elementary excitations. 
Some of the experiments are: 

• A phase separation between two paramagnetic phases - monopole liquid and 
monopole gas - as seen in kagome ice [53] 

• Magnetization dynamics (see below), whereby the time derivative of magnetiza- 
tion maps to a monopole current [54] 

• Measuring the monopole charge by an analog of the Wien effect [55]: that is, 
how the density of monopoles increases in a magnetic field which reduces the 
"ionization energy" , manifested in a reduction of the spin autocorrelation time as 
probed by muon spin relaxation [56]. 

• Thermodynamics. ..[33, 56, 57, 58] 

It should be observed, that many of the experiments basically probe the existence 
of thermally excited pseudo-charge defects as could be found in any Coulomb phase 
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system; only a few of the experiments test the special "monopole" property of the 
defects in the spin ice case (namely their literal magnetostatic interaction). 

4-2.1. Correlation length experiment A recent experiment [57] on spin ice at high 
magnetic fields was interpreted in terms of monopoles. The polarization is near 
saturation, so the only freedom is in dilute line-like excitations, of overall density 
proportional 5P = Pmax — P- These are none other than the "world lines" of Sec. 7.4 
and also are genuine "Dirac strings"; if the flux constraint is absolute, they never 
terminate. But in the presence of a small density n of defects, each world line connects 
a defect and an antidefect (and all defects are endpoints), the typical length being P/n. 
[The connccting-up is well-defined only to the extent that world lines don't touch, 
which is true if (AP)^/n is sufficiently small.] The numerical agreement between the 
diffraction width and simulation results is indirect evidence for the finiteness of the 
Dirac strings, and hence the presence of separated defects. 

5. Dynamics 

We could also make predictions for the dynamic fluctuations and response of the 
polarization field. 

5.1. Accessibility of states under updates 

As was pointed out (in Sec. 1.2), the natural update for a model with a flux is a loop 
update. [An example of such an update is shown later in Fig. 4(a).] However, with 
luck, one can run a simulation considering just minimal updates, i.e. involving the 
shortest loops (one plaqucttc, in the example of Fig. 4(a), or a hexagon in the case 
that the parent lattice is a honeycomb or diamond lattice). 

Now, in order to equilibrate a simulation, an update move must allow one to 
access every state with a given global flux, from every other state with the same flux, 
given enough steps. Minimal updates allow access in two-dimensional dimer models, 
but seem not to in the diamond-lattice dimer covering, so that non-local updates are 
needed. One method for this is a generalized "worm" update [61] (originally applied 
just to accelerate the dynamics in models for which minimal updates allow access). 
One creates a defect pair, allows the defects to diffuse till they re-annihilate, and then 
accepts or rejects the final configuration according to the usual Metropolis fashion. 

This is the typical mechanism for updating in a real system, too. That is, even 
rather dilute thermally excited defects are essential "lubricants" facilitating relaxation. 

5. 2. Dynamics of defects 

Topological defects, even when very dilute, may be crucial in the dynamic behavior for 
two reasons. First, the inherent dynamics of the Coulomb phase is presumably local, 
but a finite update may be insufficient to access all states in our ensemble. After all, 
as 1 have mentioned, the natural update follows a loop, and - given that the Coulomb 
phase has power-law correlations - one can't rule out that the updating loops have 
a power-law distribution with a tail of long loops. But if a pair of defects is created, 
random-walk for a long time, and re-annihilate, this accomplishes the same effect as 
a long update loop. (If the trails of the defects never crossed - a big "if" - that loop 
would just be the joining of those trails.) The defect pair might be created by thermal 
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activation in a realistic dynamics, or artificially as a recipe to accelerate some Monte 
Carlo dynamics. 

Defects matter for a second reason in cases, e.g. spin ice or water ice, where the 
model polarization is a physical polarization (magnetization or electric polarization) 
that couples to external fields: they control the relaxation of the polarization. Let's 
consider spin ice for concreteness, and let T > 0, so there is some density n(T) of 
thermally excited monopoles. 

Whenever a monopole (defect) moves by Ar, it changes the system's total 
polarization by APtot = QAr. Firstly, this implies the relaxation rate of the total 
polarization is the monopole current, 

Secondly, it implies the external (real) magnetic field B applies a force 

Fb = QB (5.2) 

on each monopole defect. So assuming a linear response with a mobility coefficient 
Mdrift, each defect has a drift velocity /idriftQB. Finally, the system's total 
magnetization is fjPtot, so putting it all together, the magnetization density behaves 
as 

= n(r)/idnftQB. (5.3) 

(A stochastic term is omitted). This apparently explains [54] the temperature 
dependences of the relaxation times observed years ago by [60] . 

A slab of ideal spin ice, placed in a magnetic field, responds (at least in 
simulations) just like a slab of conductor placed in an electric field. After an initial 
current, surface charge layers form which nullify the field in the interior, and the 
current comes to a halt [54]. 

The non-equilibrium recombination dynamics of monopole defects with 
antidefects after a sudden quench is worked out in Ref. [59], who suggest its use 
in real experiments. 

5.3. Wavevector (in) dependent dynamics 

Imagine first that the dynamics takes us between states that all exactly obey the 
zero-divergence constraint on P(r). Then the dynamics of P is diffusive and Standard 
dynamics would say 

^^ = r(k)4^ + C(k,i) (5.4) 

= r(k)i^Pi(k)+C(k,i). (5.5) 

I wrote a kinetic coefficient "r(k) to allow for different versions of the relaxational 
dynamics; C(k, t) is a Gaussian noise source with correlation function 2Tr(k), which 
satisfies detailed balance, and thus ensures that the steady state is the Boltzmann 
distribution implied by (3.5). Hence we get a relaxation time r,^^(k) = r{'k)K. 

The naive expectation [3] is that, since the dynamics conserves polarization, 
Tq^ oc r(k) DC |k|^, vanishing at small k. (That is the location of the static pinch- 
point singularities; since the actual variables measured are generally related to P by 
a staggering, they get offset to other wavevectors K.) However, for a vector spin 
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model, explicit analytic calculation (with a spin-conserving dynamic extension of the 
largc-n approximation and simulations (with realistic precessional dynamics) gave a 
k-independent relaxation [62], i.e. apparently r(k) ~ constant. 

This is not the full story; at T > 0, the flirx constraint is not fully satisfied, 
so there is a finite correlation length k"^ (exponential in 1/T, for a discrete model, 
or k,(TY (X T in a vector spin model), as mentioned above [see Eq. (4.2)]. Ref. [62] 
found a second contribution to the relaxation rate, tt{T), which represents a current of 
thermally excited monopoles (or the analogous excitations in a vector spin model [62]). 
The combined formula for the polarization correlations is [62] 

(P^(k,i)P.(-k,0)) oc (V - (5-6) 

+ ^2flt'' „, e[V*«^°('-)+V^-l*. (5.7) 
|k|^(|k|^ + K^) 



6. Quantum models: a quest for "artificial light" 

It is natural to ask, can we have, not just an emergent magnetostatics or electrostatics, 
but an emergent electrodynamics with separate "electric" and "magnetic" fields 
satisfying Maxwell's equations? Indeed, this can emerge in quantum versions of the 
models we discussed [11, 63, 64]. 



6.1. Set-up of quantum models 

Typically, the Hilbcrt space consists of discrete configurations satisfying the same 
constraints we discussed before, allowing the definition of a "magnetic" field via the 
same coarse-graining developed in Sec. 3. Furthermore, the model has "flip" terms of 
amplitude t^^^^ which hop you from one basis state to another. These are the same 
kind of "flips" (constrained by the flux conservation) we found in the classical (Monte 
Carlo, or physical) dynamics (see Sec. 1.2 the flip must rearrange a closed loop. 

In a quantum model it is particularly desirable the "flip" should look as simple 
as possible (at least, should be finite!) when written out in terms of spin or 
creation/annihilation operators. (Among other things, only the simpler cases are likely 
to be realized in experimental systems.) A "simple" flip means a "minimal update" 
which (in the quantum model) amounts some variety of "ring exchange". That is, it 
switches the occupancies, or reverses the spins or arrows, along the shortest loops in 
the lattice. (These are six sites in the cases of the kagome or pyrochlore site lattices). 

A natural route to the needed flip terms is to let the Hilbert space admit local 
violations of the flux constraint, but to make the Hamiltonian strongly penalize them 
by a potential energy cost of order V . The Hamiltonian also includes simpler flips of 
amplitude tpair involving two sites (spin exchange or boson hopping). Then a good 
effective Hamiltonian (low energy theory) restricts the Hilbert space to the subset of 
configurations satisfying the flux constraint, but includes loop flip terms generated 
from perturbation theory [11, 64, 68]: e.g., 

tring (X tpairV^^^ (6-1) 

for loops of length six). 
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6.2. Theory of the "Maxwell" phase 

This theory, written in terms of discrete variables, is some kind of "compact U{1) 
gauge theory". Here U{1) simply means there is one sort of "electric" charge. It is the 
electric flux which has the defining constraint I mentioned in Sec. 1, which implies a 
gauge freedom in the phases of the electric field operators. * Also, "compact" means 
the flipping term looks like cos(<i>B) where $b plays the role of the m,agnetic flux, 
and thus the magnetic flux is only defined modulo 27r instead of taking arbitrarily 
large values. Now, let's assume the $b fluctuations are rather small, jj Then events 
("instaiitons") in the quantmn evolution in which $3 wraps around its period get 
suppressed, maybe so much that we may ignore them, in which case can imagine $b 
being non-compact (taking unbounded values) as in the familiar electrodynamics. 

The next step is usually a leap to a continuum theory, by substituting for the 
"magnetic" (flip) term the simplest analytic (i.e. quadratic) term with the correct 
symmetries, with a coefficient Kb (analogous to the magnetic permeability); also, a 
quadratic "electric" term is assumed as in the classical model. If all that is justified, 
it follows that we get emergent "artificial light" elementary excitations [64, 11] with 
dispersion 

na;(q) = nc|q| (6.2) 

where he oc ^/KB/KE is the analog of the speed of light. Although linear dispersing 

excitations arc common enough in solid state (magnons in antiferromagnets, or 
phonons), these "emergent photons" are novel in that they are not Goldstone modes. 
I will call such a phase a "Maxwell" phase. 

The structure factors of such a model do have singularities at the same pinch- 
points as the classical models, but the correlation behavior has different exponents. A 
short cut to the behavior (as in any other model that reduces to harmonic oscillators) 
is the virial theorem, which says the oscillators "potential" and "kinetic" energy 
expectations each get half the zero-point energy. Thus, 

(:^|E(q)P) = 2(;..;(q)/4) (6.3) 

(the 2 is for two polarizations). The equal-time mode expectations are thus (in place 
of the classical result (3.9)) 

Se{cO ^ (E(q)„*E(q)^) « ^^^= ('^a/S - kak/3) (6.4) 

and correspondingly the correlations should decay as l/r'* in the quantum case. 
(The structure factor as a function of (q, w) is similar but proportional [64] to 
|qp/(w^ -I- c^lqp) in place of the |q| factor in (6.4).) 

I somewhat mistrust the leap to the continuum theory, if justified only by the 
symmetries of the "flip" operator. Consider the following: given a classical model 
with a classical flip dynamics, one can always devise a quantum Hamiltonian whose 
ground state has the same probability distribution as the (Boltzmann) weights of 
the classical model, via the Rokhsar-Kivelson (RK) construction [65, 66, 67]. But 

* There is not a uniform convention inthe literature whether the variables are electric and their duals 
magnetic, or vice versa. 

tt Since E and B operators are conjugate, like position and momentum of a harmonic oscillator, 
the uncertainty relation tells us the "magnetic" fluctuations can be small only when the "electric" 
fluctuations are large, as is the case in the sort of correlated-liquid/cooperative-paramagnet phase 
we've had in mind throughout this paper. 
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usually the classical system generally has relaxational dynamics, correspondingly [66] 
the quantum system has elementary excitations with lo ^ |qp dispersion; alternatively, 
if the classical relaxation time is constant (as in the model summarized in Sec. 5.3), 
the quantum excitations are gapped. 

A second worry is that, as already noted (Sec. 1.2), a rule that flips only the 
smallest loops is quite possibly insufficient to access the whole Hilbert space. 

6.3. Simulation tests for emergent electrodynamics 

Two groups have verified the emergent "Maxwell phase" [68, 69]. Their models are 
dimer coverings of the diamond lattice (in the case of [68] , an equivalent boson model, 
in which the fiip term is generated as in (6.1).) The "smoking gun" for this state 
according to [68] is correlations linear in q like (6.4), while for [69] it is an energy 
depending as )E)^ on the uniform flux density. Unfortunately, with quantum Monte 
Carlo it would be difficult to measure the corresponding dependences for the conjugate 
variable B, which ought to behave the same as E. 

7. Phase transitions out of Coulomb phases 

Now that we have (in the polarization field) something analogous to an order 
parameter, we are equipped at least to classify transitions out of a Coulomb phase 
into, say a long-range ordered phase. We might even be able to set up a field theory as 
the basis for some sort of renormalization group. But the constraints of the Coulomb 
phase necessarily mediate long-range interactions; as is long known [70] these can 
change the universality class of a transition, modifying the critical behavior to that 
characteristic of a higher spatial dimension. 

Since the long-range ordered pattern is a particular configuration, it has a net 
polarization Pq, which can be used to classify the cases. 

Zero-polarization: Pq = is zero; the transition is a symmetry breaking into one 
of several symmetry-related states, all with Pg = 0? 

Intermediate polarization: Pq 7^ so there are several symmetry-related 
directions for Pq, but each corresponds to just one kind of ordered state. 

Saturated polarizaMon: Po 7^ 0, with a symmetry breaking as in the previous 
case, but with Pq at the edge of the range of possible polarizations, so that there 
are no fluctuations whatsoever in the ordered state (so long as the constraint still 
holds) 

Furthermore, such transitions might be driven either by interactions, or by 
external fields. 

7.1. Dimer crystal transitions: zero polarization 

The most elementary case is simulated dimer coverings of a cubic lattice, with an 
interaction that favors dimers on opposite sides of a square plaquette (in d = 2, on a 
square lattice, that would induce the same columnar pattern found in quantum dimer 
models.) The ground state of this Hamiltonian is a periodic pattern having P = 
and with a six-fold degeneracy. (Three directions the dimers could line up in, and two 
choices of which layer to have dimers.) 
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Simulations of this model by Alet et a/[71] showed a continuous transition 
(whereas Landau theory would predict a first-order transition, or an intermediate 
disordered phase). Furthermore, an emergent 5*0(3) symmetry was seen in 
polarization space (by plotting the distribution of polarizations), but only right at the 
transition [72] . These properties are suggestive of the recent paradigm of "dcconfincd" 
criticality [73, 74] whereby some phase transitions that must be first order according to 
Landau's picture, instead are continuous, and the critical state possesses an emergent 
symmetry lacking in either of the adjacent phases. A flurry (still unresolved) ensued 
of theories [75, 76, 77, 78] and of foUowup simulations [78, 79] 

Before these simulations, [80] had pondered the transition from the Coulomb 
phase to an ordered dimer pattern with zero polarization, concluding it is equivalent 
to the ordering of a superconductor when electromagnetic field fluctuations are taken 
into account, which is supposed to be weakly first order. 

Ref. [78] simulated less isotropic versions of the interaction Hamiltonian. The 
basic picture of the transition is that, coming from the high-temperature (Coulomb 
gas) phase, all of these cases belong to the same universality class, an inverted XY 
transition. Those free energy terms anisotropic with respect to the polarization 
direction are irrelevant at the critical fixed point, hence the critical state has the 
SO (3) symmetry, but they are relevant at the attracting fixed point characterizing the 
ordered phase. Thus the transitions are less exotic than in the isotropic case, but in one 
case a surprising phase diagram was found: as temperature is raised, the dimer crystal 
melted into a paramagnetic phase, then at a higher temperature entered a Coulomb 
phase - re-entrant, in the sense that it has longer range correlations. (Notice that even 
if flux constraints are not enforced in the microstates, they may appear emergently.) 

Ref. [79] elaborated the basic dimer model in a different way by including further 
neighbor interactions, and found a multicritical point (i.e. the transition switches from 
continuous to first order) not so far away from the basic model simulated by [71]. 

Chalker and Powell [76, 77] applied the world-line mapping (see Sec. 7.4) to this 
problem, and massaged the result into a field theory with two complex fields. The 
symmetries related to the specific lattice imply that the terms in this field theory 
are either SU{2) invariant, or else of order 8 (hence possibly irrelevant). This makes 
plausible the observed 5*0(3) (cquivalcntly SU{2)) symmetry at the transition, which 
goes with the "noncompact CP^" (NCCP^) universality class. 

7.2. Intermediate polarization transitions 

Whenever Pq 7^ (and is sufficient to label the symmetry-broken states), it should be 
possible to describe the transition just in terms of P(r) and its derivatives. If Fq were 
not saturated, the prescription of Landau theory should lead us more or less in the 
right direction: the free energy density /(P) retains the full symmetry of the lattice, 
is analytic in P, and has global minima at the Pq values. Furthermore, rather than 
compute the exact /(P), we can usually get away with Taylor expanding it to the 
lowest order that has minima in the right places. For example, for ordering into a 
state with Pq in one of the six (100) directions, we might have 

/(P) = -ia|P|2 + - v{Pl + + P^) + ^71 VP|2 (7.1) 

When a,(3,v"f are all positive, the first two terms in (7.1) define a minimum along 
the entire sphere jPj^ = a/P; around that sphere, the third term is minimized when 
P points in any of the (100) directions. The gradient term puts a cost on the domain 
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wall between ordered domains that have chosen different ones of the six Pq directions. 
So far, this looks exactly like the Landau-Ginzburg- Wilson "action" nscd to formulate 
the e-expansion renormalization group for a 3-component ferromagnetic spin model 
with cubic anisotropy [81]. However, we are not done. If the flux constraint continues 
to hold through this transition, wc must enforce V-P = in the continuum treatment. 
We saw in Sec. 3 that this constraint generates effective dipolar interactions betwen 
polarization fluctuations in different places, and such long-range interactions were 
shown long ago to change the critical behavior of spin models. 

The case of nonzero, but unsaturated, polarizations has the fewest exotic features. 
It is realized by a Hamiltonian in which some bond directions are slightly strengthened, 
thereby selecting out a unique ground state (modulo spin rotation symmetry), so the 
system undergoes a Neel transition. As might be expected, it behaves like the ordering 
of a vector spin system with power-law interactions. In the vector-spin cases, this is 
slightly different from the related models studied in the 1970s [70], as the polarization 
field carries both spatial and spin component indices: the interaction is anisotropic 
(dipole-like) with respect to the former, but retains the rotation symmetry of the 
latter. 

7.3. Saturated polarization case 

This case has the further complication that, by the rules of our constraint, no 
fluctuations are possible at all on the ordered side of the transition. Wc thus get the 
asymmetric behavior called a "Kasteleyn transition" [82, 84, 83]. Approaching it from 
the ordered side, there are no critical divergences to signal the impending transition: 
in that sense, it looks like a first order transition. On the other hand, coming from 
the Coulomb phase side, fluctuations are allowed and in fact have critical divergences 
as one would expect from a continuous transition. Finally, states near to the ordered, 
saturated state can be expressed relative to it as a set of dilute string excitations (see 
Sec. 7.4) and that typically implies a relatively tractable behavior. 

7.4. Mapping to d — 1 quantum problem: World line approach 

Here I summarize a second analytic technique, for Coulomb phases (beyond plain 

continuum theory of the polarization field): mapping a configuration in d + 1 
dimensions to the world-lines of a set of particles in d dimensions. Thus, one is 
taking advantage of the usual correspondence between the partition function of a 
{d + l)-dimensional problem in classical stat mech and the path integral of a d- 
dimensional quantum many-body problem. The conservation law of these quantum 
particles corresponds to the defining divergence condition of the Coulomb phase model. 

The world-line mapping is ideal in the case of a transitions to a phase with 
saturated polarization, [83], since that is where the density of world-lines goes to 
zero, but was also applied to a zero-polarization case in Refs. [76, 77]. It has the 
disadvantage (compared to the polarization field theory) that it necessarily breaks the 
lattice symmetry by selecting one direction to be the (imaginary) time axis, along 
which the world-lines run. The set-up is detailed in Appendix Appendix A 

The Coulomb phase is the disordered phase of the world lines, which means that 
to every configuration found in the ensemble, you can find another configuration in 
which all the lines come out in the same place, except that two of them have got 
switched. Exactly as in the "entangled vortex liquid" phase, this signifies that the 
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bosons in the 2 + 1 dimensional quantum model are in a Bose-condensed phase - 
i.e., an ordered (!) phase with a continuous symmetry-breaking labeled by one angle 
variable. Therefore it has gapless Goldstone modes, the correlations of which decay 
as power laws in space and time. Furthermore, after one works through the mapping 
of dimer occupations to fluxes and of those, in turn, to the boson density and current, 
those power laws are precisely the dipolar correlations expected for the fluxes. 

8. Disorder 

Coulomb phase ideas can be applied to disordered cases - provided that the disorder 
does not destroy the Coulomb phase nature of the state. If defects arc very dilute, 
the basic phenomenology is that of the unchanged Coulomb phase between them, as 
discussed in earlier sections - in particular, the nonlocal spin response in Sec. 3.2.2. 
Similarly, disorder might create traps for pseudocharge defects, (Sec. 4), affecting their 
phenomenology much as impurity sites affect carriers in a compensated semiconductor. 
Here I address examples of disorder which determine the collective state of the system. 

8.1. Bond disorder 

One way to gently introduce disorder to a Coulomb phase is to weakly modulate 
the interactions in a nearest-neighbor antiferromagnet realization (say the pyrochlore 
lattice) , so the bond strength has variance [85] . Standard replica techniques gave 
a mean-field spin glass transition temperature [85] 

T^^ = ^A[x,G^y^' (8.1) 

where A4 = 6 is a coefficient derived from the adjacency matrix, and G"^ = Y- G^j, 
where is the correlation between spins i and j. Where Coulomb-phase notions 

enter is just in the evaluation of G^. This has (thanks to our assumption of small 
A) approximately the same numerical value as the pure case, for which we know the 
asymptotic behavior (Sec. 3.2), and indeed (within the large- A'' approximation) the 
exact function [40]. 

Alternatively, one considers strong random modulation but in dilute places. 
This produces local pseudospin degrees of freedom, which have effective long-range 
interactions mediated by the polarization field, [85] as we knew from Sec. 3 and 
Sec. 4. f f Thus the problem maps to a dipolar magnet dilutely occupied by spins, also 
giving a spin glass phase, as known for several decades. 

8.2. Depleted antiferromagnets 

The "hyperkagome" antiferromagnetic lattice [87, 88] is a pyrochlore lattice with 1/4 of 
the magnetic sites removed, making very regular pattern of corner-sharing triangles. 
What if we dilute at random, but place the non-magnetic substitutions such that 
exactly one site gets removed from every tetrahedron [89] . This ensemble is equivalent 
to the dimer coverings of the diamond lattice (where each dimer corresponds to a 
removed site) , and thus is highly plausible as an example of a lattice-gas realization (see 
Sec. 2). The "flux" can still be defined, and is still conserved, on such "randomly 

ft This is reminscent of Villain's pseudospins induced at wrong bonds, which had pseudo-dipolar 
interactions mediated by the Goldstone mode of an ordered antiferromagnet [86]. 
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Figure 4. Loops in ice model and constrained disorder, (a). An ice-arrow 
configuration; outlined is an update loop (reversing all arrows gives another valid 
configuration), (b). The same configuration; highlighted (pink arrows) is an 
alternating loop (ice-arrows point in alternating directions) (c). A lattice gas of 
A and B atoms (red and green circles) placed so that every vertex is surrounded by 
exactly two of each species, and represented by the same ice-arrow configuration 
as (a). (Each red atom maps to an arrow from an even to the odd vertex, and 
oppositely for each green atom.) Loops connecting atoms of the same species form 
alternating loop. Banks [17] showed that for some plausible species-dependent 
couplings (in the three-dimensional analog), the spin directions alternate along 
each loop as shown. 



depleted lattices". Hence, they still have Coulomb phases despite the disorder [89]. 
The disorder acts on the coarse-grained polarization field like a random anisotropy 
term with long-range dipolar- like correlations, since the randomness itself is a Coulomb 
phase ensemble [89]. 

8.3. Constrained disorder models: loop-based dilution 

Besides the depleted lattices, there are other examples of the (nearly unstudied) class 
of correlated-disorder ensembles that I'll call "constrained disorder", meaning the 
frozen configurations come from some "Coulomb phase" ensemble. For example, it 
was proposed by [17] that in CsNiCrFg [16], A (=Cr) and B (=Ni) spins populate 
the pyrochlore such that every tetrahedron has two of each. As noted in Sec. 2, this 
realizes an ice model. 

The spin-spin exchange constants are Jaa, Jab, aud J_bb, depending on the 
species at either end of a bond. In a realistic range of couplings, the A spins in each 
tetrahedron point along -f-n and — n while the B spins point along an independent 
direction -|-n' and — n'. Hence, the lattice breaks up into disjoint loops of A or B spins; 
in each loop spin directions alternate, and each loop chooses a staggered direction 
independent of the others [Fig. 4(c)] Thus the disorder-averaged spin correlation 
is simply the connectedness correlation function - the probability that sites i and j 
are on the same loop. In simulations [17], Cij oc 1/rfj with a = 1. 

Consider flux-loops, as mentioned in Sec. 1.2, and also the basis of world-lines 
in Sec. 7.4: flux runs in the same sense along the loop or line. One can make a 
one-to-one correspondence between configurations which have a flux-loop containing 
both sites r and r', to conflgurations with a charge defect at r and an antidefect 
at r'. (Namely, one maps the former to the latter by reversing the flux along half 
the loop between r and r'. On the other hand, if one has a charge defect, it must 
be the endpoint of two such flux lines; and if the system has only one defect pair, 
the other endpoints must be at the second defect, so the map is invertible.) This 
idea - originating with [90] - was applied in two-dimensional height models to get 
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the connectedness correlation function [47, 46, 91], and the same idea will work in 
Coulomb phases. However, the loops in question here are of a different kind: they 
run alternately with and against the flux, they cannot be reversed, so the defect-pair 
trick is not applicable. It will be interesting to see how a simple power-law can emerge 
for the connectedness correlation function. The same kind of loops are formed by 
fermions in the wavefunctions envisaged by Fulde et al [19, 20]. 

9. Conclusion 

When local constraints imply the existence of a conserved "flux" in some statistical 
model, this flux can be coarse-grained into a polarization field which is a comprehensive 
framework - not yet exhausted - to map out the model's behaviors: long-range 
correlations, charge-like defects, dynamics, and disorder eff'ects. 

In the course of the review, I've outlined the associated calculational tools: 

(1) Continuum theories, e.g. to model the analytic form of correlations in classical 

or quantum models; 

(2) The large-n limit of n-component spin models, to evaluate correlations in 
particular lattices. 

(3) "World lines" (or equivalcntly "Dirac strings") near the limit of saturated 
polarization, e.g. the behavior of "spin ice" magnets in large fields of the 
appropriate orientations. 
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Appendix A. Set-up of world-line mapping 

Here I present details of the world-line approach, applied to a variety of questions about 
Coulomb phases by Clialkca' and c;ollaborators. THe key idea is to map a c;onfiguration 
in d 4- 1 dimensions to the world-lines of a set of particles in d dimensions. Thus, one 
is taking advantage of the usual correspondence between the partition function of a 
{d -\- l)-dimensional problem 

For pedagogical purposes, my illustrations are from the 1-1-1 dimensional cases, 
although (as noted already in Sec. 3.3) these are "height" models and thus have much 
more structure than the 2 + 1 dimensional Coulomb phase models. Specifically in the 
world- line mapping, additional features in 1 -|- 1 dimension are (i). if the model is 
exactly soluble, this is how one sets it up (e.g. for the Bethe ansatz) (ii). usually the 
world lines are not crossing, so they represent hard core bosons; in one dimension, 
one can re-interpret the particles (via a Jordan- Wigner transformation) as fermions, 
and the statistics naturally implements the hard core constraint. In 2 -|- 1 dimensions, 
one is stuck with hardcore (thus, interacting) bosons. Since one is mainly after long- 
wavelength properties, the repulsive interactions are not a big concern: one coarse- 




Figure Al. Mapping configurations of a flux model to world lines. Parts (a,b,c) 
illustrate the dimer-diflfcrence construction, (a) A typical dimer covering on the 
honeycomb lattice, with arrows giving the fluxes; single arrows are uncovered 
edges, (b) A reference dimer covering, having the maximum possible flux (in the 
= y direction), (c). The difference between the fluxes from (a) and (b) forms an 
array of non-intersecting lines that always move from top to bottom, and may thus 
be interpreted as world lines of a set of particles. Parts (d,e,f) show the uniform- 
flux construction on the square lattice: again, (d) is a typical dimer covering, 
(e) is the reference flux pattern [but unlike (b), it does not represent any dimer 
covering], and (f) is the difference, showing world lines. When we slice through 
every fourth layer of bonds [dotted lines in (f )] , the lines are always paired in that 
layer, and we interpret every pair as a particle [filled circle] in a chain of sites. 
(In (c) and (f), all world lines are equivalent, but those which go with different 
particles are alternately colored red and green.) 



grains beyond the scale of the boson spacing, and describes the boson fluid using just 
(quantum) hydrodynamics. 

Our mapping must ensure that the world lines don't turn back on themselves (so 
that the particles are conserved and don't annhilate). This is done by letting the world 
lines be the difference between the actual configuration and a reference configuration 
having saturated polarization in some reference direction. If the original configuration 
had a nonzero mean flux along the reference direction, the resulting array of world lines 
has a different density from the reference; while if the original flux was transverse, the 
resulting array becomes tilted, representing a nonzero current in the lower dimensional 
quantum system. 

Figure Al illustrates two classes of this construction. Fig. Al(a,b,c) is the "dimer- 
difference" construction. Quite generally, if you take two dimer coverings of the same 
lattice and fill in edges wherever they differ, you get a set of unbroken lines or loops. If 
one of the coverings had saturated flux in a certain direction, these lines are directed 
along that same direction. This construction was used in [92] . 

I call the other construction a "uniform-flux construction" ; this is implicit in the 
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set-up used by Powell and Chalker [76, 77]. They define it for dimers on a simple- 
cubic lattice, which is oriented with the reference being a (111) axis; Fig. Al(d,c,f) 
shows the square lattice analog. The diagonal orientation of the world lines has the 
virtue that coordinate-axis directions remain symmetry-equivalent; (Of course, one 
cannot manage that on other lattices, e.g. the honeycomb; and if you wished the 
(11) type axes to remain equivalent, you would pick the other orientation.) In the 
uniform-flux construction, the flux lines are not literally world-lines, since they can 
turn momcncntarily backwards in the odd layers; furthermore each "particle" has two 
flux lines which split apart. However, if one looks at the particles only in every fourth 
layer (dashed slices in Fig. Al(d), the lines always move forwards in time (downwards 
in the figure) and always rejoin in the same pairs. In the simple cubic lattice, the slices 
should be taken every six layers, and each particle splits into three flux lines, which 
won't always rejoin but usually do. Rather than compute the exact transfer matrix 
eigenvalue \^{q), Powell and Chalker just use a function with the same behavior near 
its maximum, since they are only interested in long-wavelength behaviors (and will 
get those by mapping the d = 2 + l problem to a field theory). 
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